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I . INTRODUCTION 


The present work is an attempt to find an analytic approach to 
the oblique shock wave - laminar boundary layer interaction problem. 

It is well known that when a shock wave impinges upon a stationary 
object in the flow, there results a complex interaction between the 
shock wave and the boundary layer on the object, in which the press- 
ure rise can feed upstream along the object and the resulting inter- 
action can, for certain flow conditions, lead to separation (see 
Fig. 1), The thickening of the boundary layer caused by this up- 
stream influence disturbs the external flow and it is therefore 
impossible to compute the boundary layer independently of the 
changes in the outer flow caused by the shock. Some free inter- 
action conditions must be established to link the two regions. 

The present work attempts to find a solution in the case of a 
two-dimensional flow along a flat plate at zero incidence. 

The earliest works in this field are experimental. These include 
the works of Ackeret, Feldmann and Rott [1], Liepmann, Dhawan and 
Roshko [2], Chapman, Kuehn and Larson [3] and Hakkinen, Greber, Trill- 
ing and Abarbanel [4] . Later investigators used integral boundary 
layer solutions coupled with an appropriate pressure law. These works 
are based on the theory of Crocco and Lees {5,6]. These methods did 
produce solutions that separated, reattached and met a downstream 
boundary condition of zero pressure gradient flow, but were limited 
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to a small range of Mach numbers, Reynolds numbers and wall temper- 
atures. 

Reyhner and Flugge-Lotz [7] carried the approach a step further 
by coupling a finite difference boundary layer solution with a 
Prandtl-Meyer pressure law and successfully obtained solutions 
that included a separated region and met the downstream boundary 
conditions dictated by a zero pressure gradient flow. 

All the previous solutions rely on two basic factors: 

1. Garvine [8] has shown that when a boundary layer is 
coupled to the inviscid stream, the problem assumes an elliptic 
character in that the upstream boundary conditions cannot be com- 
pletely specified but rather that the correct solution is dependent on 
satisfying a downstream boundary condition. 

2. Interacting boundary layer solutions exhibit a nodal 
behavior in that they tend to branch away from weak interaction 
solutions towards either separation or an expansive flow (see Fig. 2). 

The advantage of the present work over the previous ones is in 
the improved treatment of reverse flow regions: the customary device 

of eliminating or modifying the time-like terms in the streamwise- 
momentum and energy equations in this region [7, 9, 13] has been 
eliminated by adding a time-like term to both sides of these 
equations in the backflow region, thereby avoiding the change of 
the basic nature of these equations, as called for by the previous 


methods . 
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The main advantage of using the method employed in the present 
work is that it provides a solution to a given case in a short time 
(usually about one minute of computer time using the UNIVAC 1110) as 
opposed to considerably longer time when a Navier-Stokes procedure 
is employed using present technology computers. This has a strong 
bearing on cost considerations, yet the accuracy of a solution as 
produced by the present method is good enough for engineering pur- 
poses . 

The present effort assumes that the Prandtl-Meyer relation 
adequately represents the external flow and the free interaction 
condition between the boundary layer and this external flow is 
through coupling the two regions along the displacement thickness 
surface by a relation between the external flow turning angle and 
the stream-wise derivative of the displacement thickness. Implicit 
in this work is the assumption that the boundary layer equations are 
applicable throughout the computational regime, which excludes large 
separation bubbles, where the normal pressure gradient is no longer 
negligible. 

There are three par- aeters whjch are of fundamental importance 
in this problem: the starting point of the interaction, the location 

of the idealized shock, and its strength. Knowledge of the latter 
two of these quantities does not permit the knowledge of the third, 
and hence iterative calculations must be performed to find the third 
parameter, namely, the starting point of the interaction. The 
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determination of this parameter is done through meeting the dovmstream 
boundary conditions, which are: (a) weak interaction far enough down- 

stream of the shock impingement location; (b) a downstream pressure 
which is compatible with the one predicted from purely inviscid con- 
siderations. These two conditions necessitate the employment of two 
iteration parameters, such that only a unique combination of them 
yields the correct downstream conditions. Adopting a combination of 
the work of Reyhner and Fluggvi-Lotz [7] and of Dwoyer [9], the two 
iteration parameters employed in the present work are the starting 
point of the interaction (assuming knowledge of the shock impingement 
location) and the second derivative of the displacement thickness at 
the shock impingement point. 

The boundary layer equations are solved by a linear implicit 
finite difference technique [10]. 

Special treatment is applied for the backflow region, consisting 
9 u 

of adding a (-2pu -^) term to both sides of the streamwise momentum 

9T 

equation and a (“2c^ pu — ) term to both sides of the energy equation, 
such that diagonal dominance of the resulting tridiagonal matrix is 
assured. This method is more accurate than the one used by Reyhner 
and Flugge-Lotz [7] and later by Dwoyer [9], which consisted of re- 
placing the time-like terms in the streamwise momentum and energy 
equations by one-tenth their absolute values. No such artificial 
terms appear in the present work. The stabilizing method employed 
here is based on the successive iterations performed on each down- 
stream station, so that by the time a station is converged (i.e., 


\ 
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the free interaction condition is met) the streamwise momentum and 
energy equations are restored to their correct form to within a small 
prescribed tolerance. 



II. LIST OF SYMBOLS 


a 

C 


f 

g 

G 

JJ 

JD 

K 


K1.K2 


L 

M 

m 

P 

P 

Pr 


speed of sound, a = 

sec 

Chapman-Rubes in constant (eq. 21) 


skin friction coefficient, C, ... ^ v 

f (l/2)p u “ 

00 00 

specific heat at constant pressure. 


specific heat at constant volume, [ 


cal 

kg.'C 


1 

^kg.®C^ 

] 


exponent in power law for the viscosity 
exponent in power law for the thermal conductivity 
function of y used in Prandtl-Meyer calculation, G = 
counter ^or mesh points in the x-wise direction 
J at the last computational station, x = x downstream 
counter for mesh points in the y-wise direction 
coefficients in the x^-iteration procedure 


thermal conductivity [- 


cal 


r] ; also counter for mesh 


"m*sec* 

points in the y-*wise direction (will be clear according to 
context) 

reference length, [m] (L = x ) 

o 

Mach number 

number used to divide AX to reduce mesh size in x-direction 
r N T 

pressure, [-y] 
m 

nondimens ional pressure, P = 

Po 

pressure at the dovmstream station 


Prandtl number, Pr 


y c 
-P . P 


6 



< H ^5 p 


7 


q 


R 

R 

e 

X 

t 

T • 


n 


V 

V 

X 

X 


heat flux [ 


cal 

2 

m *sec 


] ; also quotient in power series for AX in the 


case of variable mesh size (clear from context) 


gas constant, [ 


m • N 


Reynolds number based on x, R 


P u 
o o 


temperature, [*^C] 
nondiraensional temperature, T = 

o 

temporary temperature vector used in backflow region 
calculation 

intermediate solution vector for the temperature at a 
streamwise station, improved by each iteration at that 
station 

streamwise velocity, 


nondimensional streamwxse velocity, U = — 

u 

o 

temporary streamwise velocity vector used in backflow 
region calculations 

intermediate solution vector for the streamwise velocity 
at a streamwise station, improved by each iteration at 
that station 
normal velocity, [— ^] 

p vL 

nondimensional normal velocity, V = 


streamwise coordinate, [m] 
nondimensional streamwise coordinate 


X = -- ^ , ' Z 
p u 
o o 
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X streamwise location of the beginning of strong interaction 

o 

X , streamwise location of the idealized shock impingement on 

Sn 

the plate 

y normal coordinate, [m] 

y 

Y nondimensional normal coordinate, Y = ^ 


Greek Symbols 


6 

Y 
6 

6 * 

AX 

AY 

y 

V 

T 

P 

G 

n 


pressure gradient, 3 = ^ 

specific heat ratio, y = (= 1.405 for air) 

boundary layer thickness, [m] 

boundary layer displacement thickness, [m] 

streamwise mesh size 

normal mesh size 

viscosity, — ] 
m*sec 

Prandtl-Meyer angle, [^ad] 
flow turning angle, [t^ad] 

I. r N T 

shear stress, [— 
k ™ 

density, [-^] 
m 

tolerance 

nondimensional normal coordinate used for initial profiles. 




Subs c ripts 


o 

00 

a 

c 

d jdownstr . 
f 

g 

1 

INV 


evaluated in the free stream at x 

o 

evaluated at the outer edge of the boundary layer 
adiabatic 

a computed value (i,e., 5*'*) 

evaluated at the last streamwise station 

evaluated at a far downstream or final station 

a guessed value (i.e., 6*”) 

g 

a value for an iterated quantity at the ith iteration 
a quantity evaluated through inviscid considerations 


initial 

j 

k 





w, wall 

X 





a quantity evaluated at the initial streamwise station 
evaluated at the jth streamwise station for any k 
evaluated at the kth normal station for any j 
evaluated at the (j,k) mesh point 

evaluated at the streamwise station corresponding to x 

o 

evaluated at the wall 

derivative with respect to x. Also, evaluated at x (clear 
from context) 
evaluated at x = x 


o 

evaluated at x = x , , 

shock 


(shock impingement point) 
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Superscripts 

’ total derivative with respect to the appropriate independent 

variable 

* nondimens ional quantity, used for p, k and y 

nondimens ional quantity used for 6 and 6* 



III. GOVERNING EQUATIONS 


The flow field is divided into a boundary layer and an external 
flow. The external flow is handled as a Prandtl-Meyer edge con- 
dition for the boundary layer. Although the Prandtl-Meyer relation 
is isentropic and correct only for expansions, it can be used for 
weak compressions since the change in entropy is proportional to the 
cube of the pressure change through compression waves, and this is a 
very small quantity when proceeding from one station to the next 
while marching streamwise through the mesh points. The presence of 
the shock is compensated for by changing the reference Prandtl-Meyer 
angle w en the computation procedure arrives one station downstream 
of the shock impingement point. The amount of change in the reference 
angle is automatically determined by the combined strength of the 
assumed shock and expansion fan and is corrected while the latter is 
iterated during the computational procedure. 

Figure 1 shows the physical phenomenon treated by the com- 
putational procedure, indicating the important locations and geo- 
metrical features, is the location where the compression waves 
first appear, and consequently where the pressure begins to rise from 
its upstream value. 

This compression turns the flow away from the wall with sub- 
sequent separation possible. This is followed by the combination of 
an oblique shock and an expansion fan at x . This is the location 
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where the boundary layer edge is most removed from the wall, and 
where the flow turns most sharply towards the wall. The distance 
between and is the measure of the extent of upstream influence 
present; the stronger the shock the larger this distance, and hence 
the larger the resulting backflow region. Further downstream, a 
second set of compression waves turns the flow back parallel to the 
wall, followed by reattachment and gradual reversion back to a weak 
interaction region. Since the present work does not incorporate a 
normal momentum equation (and hence assumes mild flow deviations 
from the flat plate and constant pressure normal to it) only small 
separated regions can be treated by the present approach. 

III. a Boundary Layer Equations 

The basic equations for the compressible, steady, two dimensional 
boundary layer, assuming a constant Cp, are 

A perfect gas is assumed, for which 




( 4 ) 
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Ppwer-law relations are used for both the viscosity and the 
thermal conductivity: 

where u and k are respectively the viscosity and thermal con- 
o o 

ductivity at the reference temperature t^. A value of 0.76 has been 
used for both f and g. 

The boundary conditions are 
u (x,o) = 0 

V (x,o) = 0 (for impermeable wall) 

V (x,o) = V (x) (for suction) 

constant wall temperature case 

constant wall heat flux case 

where u (x) and t (x) are furnished by the solution of the inviscid 

00 00 

outer flow, to be discussed later. 


u (x,“) 
t (x,o) 
t (x,“) 


w 

= u^(x) 
= 

= t (x) 


or 


1 

S^ly-o ' 


t (x,“) = t (x) 



14 


Ill.b Dimensionless Forms and Reference Quantities 

The dimensionless variables employed in this work are chosen so 
as to achieve two goals: (a) creating flow variables that are 

most convenient for the finite difference procedure in which they are 
to be used; (b) creating as few parameters as possible, to maximize 
computational convenience. 


With these in mind, the following dimensionless variables are 
employed: 



where all quantities with subscript o are evaluated in the free 
stream at x^, the point where the strong interaction starts, and L 
is a reference length in the direction normal to the wall. 

With these new variables, the dimensionless equations now 
become 




(2a) 




15 


(T-I)Mo/**(^ 

p^f*T 


ye*=T-# 

k*= 7 -^ 


where the free-stream Mach number at x is 

o 

M = (4a — 

met, 

and the Prandtl number at the same location is 

The corresponding dimensionless Boundary Conditions are: 


U(X,0) 

V(X,0) 


= 0 



0 (impermeable wall) 

y, (X) (suction) 

w 


'^00 

U(X,») = = U (X) 

U CO ' 


(3a) 

(4a) 

(5a) 

(6a) 


For constant wall temperature 


T(X,0) = Tw 

T(X,») — = T^(X) 

o 

while for constant heat flux q 


-k* 


Hi 

3Y ' y=0 


= _aL_ 

k t 
o o 


T(X, ■») = T (X) 

00 


Note : In all subsequent equations the asterisk used to indicate 

dimensionless density (p*) , thermal conductivity (k*) and dynamic 

viscosity (u*) will be dropped, with the understanding that these 

variables are dimensionless as defined above. 

The following quantities are given in the free stream at 

X = X : M , p , p , t and k . Based on these, the reference 
o o o o o o 

quantities are defined as follows: 


Velocity: a.= M, Jf9.fi 

Static pressure: 


The reference length, L, is chosen as the distance from the leading 
edge of the plate to the location of the start of the strong inter- 
action, namely 

L = X. 

which leads to the following dimensionless coordinates: 
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where 



III.c Edge Conditions 

The edge conditions necessary for the solution of the governing 
equations are calculated from the Prandtl-Meyer relation which is 
strictly valid only for isentropic expansions. However, the 
Prandtl-Meyer relation can be used for weak compressions since the 
change in entropy is proportional to the third power of the change in 
pressure for compression waves, and this entropy change is very small 
when proceeding from one mesh point to the next one while marching 
downstream during the computational procedure. This approach has 
been found also by previous authors to work reasonably [7,9,13]. 

The main difficulty in the solution of the governing equations (1-6), 
is that the pressure, p, is one of the unknowns, since in the case of 
a strong interaction the viscous and inviscid regions mutually 
affect each other and the pressure field is a result of this inter- 
action. This is the basic difference between strong and weak inter- 
action cases, since in the latter the pressure field is dictated by 
the inviscid flow, so it becomes a known variable as far as the 
viscous layer is concerned. In the present situation an additional 
equation is necessary which expresses the mutual determination of 
the pressure field by the viscous and inviscid regions. This is 
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done by coupling the two regions along the displacement surface (see 
Appendix Al) by the relation 

= (7) 

where the index » refers to the conditions at the outer edge of the 
boundary layer. ^ is now used in the Prandtl-Meyer relation through 
the Prandtl-Meyer angle v: 


^ = ( 8 ) 

where 


^ J tv/ (^^9^1) - (9) 

and is the reference Prandtl-Meyer angle. The value of is that 
for M^ until one station downstream of the shock impingement point, 
where it assumes a new value that depends on P^ through the combination 
of incident shock strength and expansion fan such that it also guaran- 
tees continuity of pressure across the shock impingement location. 

This new value of is updated during the iterations performed on the 
shock strength plus expansion fan, as part of the computational pro- 
cedure. Once equation (9) is solved for M^ (see Appendix B) , the edge 
conditions may be evaluated as follows: 



( 10 ) 


( 11 ) 
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( 12 ) 



and the pressure gradient 


(13) 


T 


Mo 1 





(14) 


In equation (14) the term jr appears, which has the following 

CfX 

value: 



The first term on the R.H.S. vanishes for a flat plate. It should 

be noted that the second term contains both 6* and x . Hence 

o 

changes in both influence the computation of the edge flow para- 
meters, independently of each other. This has a strong bearing on 
the method of iterations employed in the present work which is 

basically a double iteration loop where 6* and x each change while 

o 

the other is frozen in its latest value. 



IV. METHOD OF SOLUTION 


i 


li 

(■ 


IV • a Difference Equations 

The basic equations, eqns, la - 6a, will now be rewritten in 
finite difference form. A linear implicit form has been chosen, 
having the distinct advantage that all equations are linear in the 
various unknowns, thus at each step the momentum equation may be 
solved for the streamwise velocity values at the new station, 
followed by the energy equation yielding the values of the temper- 
ature at this new station, then the density is obtained from the 
equation of state, and based on these the continuity equation is 
solved to provide the normal velocity at the new station. Finally 
the viscosity and thermal conductivity are found and the procedure 
may be advanced one step downstream. Backward differencing is 
employed for the streamwise derivatives and central differencing is 
used for the normal derivatives in both the streamwise-momentum and 
the energy equations, while backward differencing is used in both 
directions in the continuity equation. The differencing for both 
the momentum and energy equations is correct to order (AY)^ in the 
normal direction and to order AX in the streamwise direction, whereas 
for the continuity equation the differencing is correct to orders 
AY and AX. Since the truncation error is of second order in the 
streamwise direction, a relatively large number of mesh points has 
been employed, resulting in a fine mesh. The difference equations 
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IV . c Matrix Form of Equations 

The momentum equation Is rewritten in the form: 
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where 


~ 2CAyj ^ 

A ^ p^Ai^ 

Ik AT ^ ^ (Ay)*’ 

r\ — -J^jr A Vj//?, ^ < V jj } 

^ca/; l^yr 

^j’t* L_ Pj»/ ~~ fif 

4x 4Z 


This matrix equation is used to solve for U - . using a tridiagonal 

J “^1 > K 

matrix inversion technique. 

The energy equation is written as: 


Ik. 

).k f (^•i 

Pm. (Ay)^ ’ P/t 

£ kj,it 
fit 

Ti*ok 


^Ayr 




LiidiV 

[ tCA 

[ AT 

iijJ-r- 

TfSyP^ - 


,e 


Uit/.k-i 

- zM 


4. 


JZi. 


J'" AZ 


(3c) 
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Wfritten out for all k for the case of constant wall temperature, we 


get 
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4 


This is solved for ^ using, again, a tridiagonal matrix 

inversion. Now we can find from the perfect gas law 


P _ — Pif/ 

ru,,k 

Next the continuity equation is solved to yield 


(4b) 






(2c) 


This solution is started at the plate (y=0) where V.,, =0 (or 

3+1 »o 

V..- = V (X)) and marched out to the free stream. Finally, the 

j+l,o w ’ 

viscosity and thermal conductivity for the (j+l)st station are 
computed : 


/ jH , k ” ('^■*■0 


(5b) 




(6b) 


This completes the solution for the (j+1) st station and the method 


is now applied to the next downstream station. 



V. COMPUTATIONAL PROCEDURE 


V.a Initial Profiles 

The initial velocity profile is generated by an iterative pro- 
cedure that is started with the Pohlhausen polynomial for the flat 
plate laminar boundary layer. 

= 2 ^- 2 ?’*^* ( 15 ) 

where 



In terms of the Y -coordinate. 



where 6 . ^ is tentatively taken at its Mach zero value 

inxtxal 

5-0 

The initial station, x. _ . ,, is the predetermined station where 

initial 

the computational procedure starts . Hence 



(15a) 


(15b) 
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with k going from zero up to corresponding to the last 

mesh point in, the normal direction. 

The initial representation of the temperature profile uses the 
Crocco relation between temperature and velocity for a flat plate 
with zero pressure gradient which is strictly valid only for Pr = 1 

( I- ( 16 ) 

which becomes, in finite difference form, 




•» 2 


Mo ( 16 a) 

Since these profiles do not correspond to the situation in the flow 
field (M 0, Pr 5 ^ 1) at our initial station, they are refined 
iteratively by using them to solve the equations of motion for one 


step downstream with zero pressure gradient, and then using the 
resulting profiles as new initial profiles. This procedure is 
repeated 10 times to ensure realistic initial profiles that are 
solutions to the difference equations and thereby allow smooth 
marching downstream from them. 


Treatment of Reversed Flow 


The matrix U-equation was found to be 
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r N 


r \ 

ft a, 


Vm. 


+. 






/j 







““ 


► 

«— .M 


— 


— 




— 


— 

•^.-1 ^., Jlw^ 

1 




Md 






where 


J/ 

2C4y; 


(17) 








4Z 




(18) 
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In order to maintain numerical stability of this tridiagonal equation, 
diagonal dominance (DD) must be preserved. This means that the 
following must be maintained: 




or 




AX 


> 0 


This shows that DD will be lost when U. . 1 0, that is, in a 

reversed flow region (see Appendix A2) . 

In order to avoid this problem, the usual method employed by 

several authors (C.F. Ref. 7) is to drop the time-like term 

pu*u^ for u _< 0, which provides marginal stability: 

\ \ = 0. or, to assure that \ \ \ > 0. the term 

pU"U is replaced by 0.l|pu|u wherever pu < 0. This means that 
P- uU. , * 

"1 K T Ic 

the term > in the expression for eq. (17), is replaced 

... 

^ • Physically this implies adding a small artificial 

convection term to the momentum equation. The claim usually made is 
that this creates only a small additional error as compared with 


that caused by the neglect of the P’^'^ term, and it guarantees DD. 
However, both of the above methods actually change the basic 
nature of the momentum equation in the backflow region, and this 
should be avoided . 

The method employed in the present work is to add 
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2p^ k 

AX ^^j+l,k sides of the momentum equation, eqn. 

(Ic), whenever ^k 0* By doing so it is guaranteed that the 

momentum equation ^mains unaltered in the backflow region. There 

is a difficulty in doing this, however, since the term added to 

the R.H.S. of the equation in its matrix form, i.e., to (J) , 

Ic 

eqn. (18), contains the term which is not known (the 

■'^®ctor is the unknown of the matrix equation). This 
problem is solved by taking advantage of the fact that there are 
iterations performed at each J-station. Thus, a term 
2p. ,U. , 

~ AX \ added to (J>k for U^^k < 0. where Uk is the solution 
'^j+l,k from the previous iteration on the J+lst station. 

This term becomes closer and closer to the correct value, U, , 

j+1 f k 

with each iteration, and by the time the last iteration on the J+lst 
Station is reached, the following exists: 


ro 

I i/k 

' 0i..k 


<S 


and hence ^jj,is within e from thereby guaranteeing that the 

momentum equation at each J-station is correctly represented whether 
in the backflow region or not. 

By employing this approach, the effective result""is that the 

k^i k P- v|U. . I 

term ' is replaced by for U^^k < 0* thereby 

guaranteeing DD, and yet no alteration of the momentum equation is 
caused as a consequence. 
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A similar situation holds for the energy equation: 

The matrix form for the energy equation was found to be 
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“ ^(ly) 


i_k 


?5l"c5yF 






Again, to maintain DD it is necessary that 


which means that the 

■Ai K ^j/h. X o 


following must hold: 


i k i k 

Hence the term ( tx ' ^^j+1 k 

the energy equation, eq. (3c), so as to keep the energy equation 
unchanged in the backflow region. This time a temporary Tj^ vector 
is needed to multiply the above term when adding it to (J>^, eq. (20). 
The solutions Tj^^ ^ of the successive iterations on each J station 
are used for this purpose, and by the time the last iteration is 

reached, station J+1 is converged and the relation 

I IV I 




<8 


is guaranteed, thereby keeping the energy equation unaltered in the 
backflow region and yet maintaining DD there. 


V.c Iteration Logic 

As can be seen from the Prandtl-Meyer equations employed in the 
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present work, eqs. 7 - 14a, the edge conditions are based on 6* and 
its 1st and 2nd derivatives • Since 6* is dependent on the velocity 
and density profiles, which are themselves dependent on the edge 
conditions, it is obvious that an iterative approach is necessary. 
Hence the computational procedure is as follows; 

If 

1. An initial guess is made for 6*, denoted by 

6 ^. 

g 

2. This guess is integrated numerically to 

obtain (See Appendix Cl) 

g 

3. 6* and 6** are used to obtain the boundary layer 

edge conditions through the Prandtl-Meyer procedure, 
eqns. 7,8,9,14 and 14a. 

4 . The edge conditions are used to integrate the 
boundary layer equations at the new x-wise station. 

5. From the new profiles, a value for 6* is 

computed, denoted by 6*. 

c 

6. 6* is twice numerically differentiated to obtain 
6*". (See Appendix C2) Since a three-point back- 
ward differentiation is used here, a numerical 
problem exists at station J +1, where the pro- 

^sh 

cedure smoothes out the jump in 6* that is induced 

J . This may account for the local instabilities 
^sh 

which occurred at station J . +1 and continued for 

sh 

a few stations downstream of it, as will be seen in 
the figures. 

7* 6**' and 6*” are compared. If they do not agree to 

g ^ 

within a required tolerance, a new guess for 6*” 

g 

is made, using a Newton—Raphson procedure (see Appendix 

C3) , and a return is made to step //2. When 6*" and 

c 

6*” agree to within the specified tolerance, the 

o 



34 


present x-wlse station is considered solved and 
the whole procedure is marched pne step down- 
stream, with the last (6*”)^ used as the first 

g J 

guess for 6*" at the J + 1st station. 

The above procedure begins with the supplied initial profiles 
and marches downstream until no further changes in the pressure field 
occur, indicating that the shock wave— boundary layer interaction is 
over. 

Two distinct cases may occur at the downstream region (i.e., 
far enough downstream for the effect of the shock to diminish); 

a. the pressure does not acquire a steady value, but, 
rather, either goes rapidly to zero to increases rapidly; 

b. the pressure has a steady value, which, however, is 
different from the one physically occurring (determined by the shock 
strength) . 

The present work employs two distinct procedures to correct 

these behaviors and to pick successively better branches. 

In case (a) , the branching solutions are controlled by fixing 

6**' at the station corresponding to the shock impingement location, 

'^x Appendix D) . The perturbation in 6*** is created by using 

sh 

the value at the last station before J , multiplied by a parameter 

^sh 

called PARDEL. This in effect means introducing a pressure gradient 
discontinuously at this station, the^reby triggering a branch. The 
successive values for 6*" (through PARDEL) are determined by a 
Newton-Raphson procedure such that the correct 6 *” leads the com- 
putational procedure to a level branch, i.e,, constant pressure at 



35 


the downstream weak interaction region (see Appendix C4). Figure 2 

shows the effect of 6*” on the branching behavior. Along with this 

sh 

branch control, v is also changed at J + 1 to keep a constant 

^sh 

pressure across J , where the shock and accompanying expansion fan 
^sh 

are located. The iteration procedure described above (steps 1-7) 

is not performed at station J to keep the discontinuity in press- 

^sh 

ure gradient introduced by FARDEL. 

In case (b) , the branching solutions are controlled by fixing 

a new value for x (see Appendix D) . The closer x is to the point 
o o 

where the ideal shock impinges on the boundary layer, (x , ), the 

sh 

weaker the resulting interaction and the weaker while the further 

is upstream of the shock impingement point, the stronger the 

resulting interaction and the stronger P^, Hence too compressive 

branches are corrected by moving x downstream, closer to x , , 

o . sh 

while too expansive branches are corrected by moving x upstream, 

further away from x^^ (see Appendix E) , Changing x^ changes Re^ , 

of course. Figure 3 shows the effect of x^ on the branching 

behavior, indicating that minute changes in x suffice to impose a 

o 

large change in P^. This behavior has been indicated also by 

Reyhner and Flugge-Lotz [7]. While an x^ iteration is performed, 

6*” is frozen at its latest value, and vice versa. Usually only a 
sh 

few 6*'* -type of iterations are necessary to obtain a steady 
’'sh 

downstream pressure. Once this is attained, all further iterations 

are of the x -type, with 6*" frozen' at its latest value. 

’'sh 
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The accompanying flow chart (Fig. 4), explains the computational 
procedure, and Fig. 5 shows how 5*" varies with x. This figure 
shows regions of both strong and weak interaction within the same 
case. For the weak interaction it is possible to use the relation 
for the displacement thickness in the case of a compressible laminar 
flow with Pr = 0.72 [20] as a check on the numerical 

results: 


=/^[l.7208 + I) 

0.287<5(!f'-OMl] 

J D 


( 21 ) 


where C is the Chapman-Rubes in constant C = — ” which is unity 


P U 

•^00 CO 


in 


our case, since T =1.0. 

w 


Inserting numerical values into equation (21) and using the 
t 

fact that — = 1, yields, for the particular Re of Fig. 5 (in 


sh 


dimensional form - all lengths in meters) : 


S'* -8.S^37*\0'*J>r 


(22) 


and 


^’=-2.1574 x/o‘V 




and, with 

X=Ct-;)4Y 

equation (23) becomes 


(23) 


(t-i) 


-5/4L 


(24) 
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where the numerical value of Ax for the particular case of Fig. 5 
has been used . 

It is observed that the strong interaction case sharply deviates 
from the weak one between J = 60 and J = 150. The oscillatory 
behavior between J = 60 and J = 85 is attributed to the inherent 
numerical instability that results from numerical differentiation. 

As may be observed in Fig. 12 the boundary layer thickness (and 
hence also 6*) exhibits a mild oscillatory behavior over the region 
corresponding roughly to 60 < J < 85. (Emphasized by heavy dots.) 
When such a curve is twice numerically differentiated, the 
result may lead to unstable oscillations, since the 
nature of differentiation is to amplify any existing instability. 

However, the iterative procedure used in the algorithm always 

managed to stabilize the computations based on 6*” since the effect 

of the iterative procedure is to produce an average trend line of 

6*" in the oscillatory region. Downstream of J , 6*” > 0 (dis- 

o 

regarding a few isolated points where 6**’ < 0 as a result of local 

numerical instabilities downstream of J , where reattachment 

sh 

occurs) , which is characteristic of strong interactions (see also 
Fig. 19), Around J = 115 the strong interaction begins to decay 
and approaches the downstream weak interaction region asymptotically. 



VI . COMPUTED RESULTS 


Several examples have been run to test the method used in the 

present work. Common to these examples are the following flow 

parameters: M = 2.0, T = 1.0, Pr = 0,70, and Re = 3.89 x 10^, 

o w sh ’ 

corresponding to an altitude of 16,500 meters. The only parameter 

that was varied from one case to the next was the overall pressure 

ratio, P^, which in effect, selected the shock strength. Examples 

have also been run with an adiabatic wall (T =1.67) but became uni- 

w 

stable for shock strengths corresponding to ^ 1.025. 

X. 

All computations were initiated at — = 0.50, with 

X u 
sh 

constant step sizes in both X and Y directions, the magnitudes of 
which changed during the execution of the program according to 
the current value of x^, such that the final step sizes were 
determined by that x^ which produced the correct branching solution. 
The mesh size consisted of 200 points in the streamwise direction 
and 72 points in the normal direction. Increasing the number of 
points above 200 in the streamwise direction did not produce any 
change in the computed results. 

Figure 6 shows the development of the streamwise velocity pro- 
files over the flat plate. The separation develops smoothly but 
the reattachment is abrupt, with the first reattached profile 
exhibiting a reduced velocity region in the upper part of the sub- 
sonic section of the boundary layer (refer to corresponding profile 
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in Fig. 8). Since this is the profile iiranediately downstream of 

shock impingement, the reason for its shape is perhaps the local 

favorable pressure gradient that occurs at this location (Fig. 13), 

which in turn may be the result of the numerical calculation of 

6*" that smoothes out the jump in ^ at J . The profiles further 

sh 

downstream do not exhibit this bulge but their shape, possessing an 
inflection point, indicates a gradual recovery from this profile. 

The fourth profile after reattachraent exhibits a small reversed 
flow near the wall. This tendency of an attached flpw to re-separate 
has been observed also by Reyhner and Flugge-Lotz [7] and attributed 


to the combination of the Prandtl-Meyer relation and the assumption 


— = 0, both of which are deficient in the reattachment zone. The 
method used here to eliminate this problem was to replace eqn. (14) 


beyond reattachment by a simple backward difference. 


1+1 1 


Figure 6a is an enlarged detail of Fig. 6, showing the separation 


region. 

Figure 7 is a comparison of streamwise velocity profiles for 
the cases of incipient separation (P^ = 1.12) and of a separated 
region (P^ = 1.31). It is observed that the thickening of the 
boundary layer is considerably more pronounced in the case of the 
separated region, as expected for a stronger interaction (stronger 
shock) . 


Figure 8 shows the development of the Mach number profiles for 
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?£ = 1.31. It is observed that the subsonic part of the boundary 
layer is roughly its inner 30% upstream of the shock impingement 
region and this changes to approximately the inner 40% of the bound- 
ary layer thickness dovmstream of this region. As the shock 
impingement region is approached, the subsonic layer increases up 
to about half of the boundary layer thickness, and decreases again 
after the shock impingement (see also Fig. 12). 

Figure 9 shows the behavior of the normal velocity profiles. 

It will be observed that a part of the strong interaction region, 
starting at J=82 may be subdivided into five regions: 

Region I : 82^j£90, with negative normal edge velocity; 

Region II : 91^J^97, with positive normal edge velocity; 

Region III : 98^J^101, with strong negative edge velocity; 

Region IV : 102^J^109, with strong positive edge velocity; 

Region V ; J^llO, with negative edge velocity, gradually 

diminishing to zero. 

The behavior in Region II is a consequence of the flow over the 
separated region, which has a component directed away from the wall. 
The behavior in Region III stems from the rather sharp velocity 
component directed towards the wall at the region of reattachment, 
while the behavior in Region IV is a result of the sharp turning 
of the flow away from the wall so as to become parallel to it 
downstream of reattachment. No explanation has been found for Regions 
I and V, except that the gradual diminishing of normal velocity in 
Region V is a demonstration of the slow switch to weak interaction, 
as can clearly be seen in Fig. 5. Since the strong interaction 
region is noted for its high irregularities (Fig. 5) and since the 
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normal velocity calculation is very sensitive to these irregularities, 
it is not helpful to analyze the normal velocity profiles in great 
detail. 

Figure 10 shows the development of the static temperature 
profiles in the case of a cooled wall. The dissipation in the 
boundary layer causes a rise in the temperature, which is then re- 
duced near the wall by the cooling (T = 1.0). 

w 

Figure 11 shows' the development of the stagnation temperature 
profiles. The reduction in stagnation temperature through the 
boundary layer is caused by the conduction of heat from the bound- 
ary layer into the wall (cooled wall case) . It is observed that 
the reduction in stagnation temperature is less in the separated 
region than in the attached regions indicating a lower heat transfer 
rate for the separated region than for the attached flows. 

Figure 12 shows the development of the stagnation pressure 
profiles. The decrease in stagnation pressure is due to reduced 
velocity in the boundary layer. Note that the decrease is sharp in 
the supersonic part of the boundary layer and becomes smaller and 
smaller through the subsonic layer, until, at the immediate vicinity 
of the wall the stagnation pressure approaches the local static 
pressure. In the separated region, because of the small recirculation 
velocities, the stagnation pressure retains nearly its wall value 
further into the flow. Also Indicated in Fig. 12 is the behavior 
of the boundary layer thickness over the plate (dotted) . The most 
upstream part behaves in the manner of a flat plate profile. The 
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boundary layer thickness then shows a moderate increase downstream 
of becoming sharper as the shock impingement location is 
approached and attaining a peak in the vicinity of shock impinge- 
ment. Downstream of reattachment there is a rather sharp decrease 
in boundary layer thickness, forming the well known "neck", and 
further downstream a flat-plate-type of profile is attained again, 
indicating transition from strong to weak interaction (see also 
Fig. 19). The heavily dotted points indicate the regions where 
the boundary layer thickness profile wiggles the most (as a result 
of numerical instabilities) and the region dpwnstream of is 
causing the oscillatory behavior of 5*" at that zone (see Fig. 5). 
For comparison, a compressible boundary layer thickness profile on 
a flat plate has also been included (dashed line) . The computed 
thicknesses are slightly below the theoretical, but close enough so 
that the computed thicknesses are credible. 

Figure 13 shows a typical static pressure profile over the 
plate, as produced by the program. It starts with a constant value 
corresponding to the upstream boundary conditions, followed by a 
relatively moderate increase starting at x^. This continues until 
a short plateau is attained downstream of the point where the shock 
impinges. At the second compression region an abrupt jump in 
pressure occurs, followed by a constant value all the way down- 
stream, corresponding to the downstream boundary conditions. The 
bulge seen immediately before the plateau is a.^result of slight 
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pressure drop just downstream of shock impingement and occured in al- 
most every case in the present work. This pressure drop is perhaps 

due to the local numerical instability introduced by 6*” at J +1, 

dP ^ 

that smoothes out the jump in — | , and may be the reason behind 

the shape of the first reattached profile (see Fig. 6). 

Figure 14 shoes a typical skin friction profile over the flat 
plate, as produced by the present program. The initial flat plate- 
type of behavior gives way, starting at x^, to a rather sharp decrease 
in leading to separation, followed subsequently by reattachment, 
and further downstream the profile again attains a flat plate-like 
shape. The wiggles in the separated region are due to local numerical 
instabilities (associated perhaps with the treatment of backflow 
regions) that plagued all cases which did not incorporate suction. 
(However, see Fig. 17.) This shape for the profile in the separated 
region has been obtained also by Reyhner and Flugge-Lotz ([7], Fig. 

15), including, particularly, the relatively large negative peak 
immediately prior to reattachraent . The large positive bulge after re^ 
attachment is a direct consequence ofthe shape of the first reattached 
profile near the wall. As may be seen in Fig. 6, this profile has a 
large slope at the wall, causing the observed overshoot in C^, This 
bulge dies out further downstream as the profile asymptotically attains 
its flat-plate value. 

Figure 15 shows the original profile for the case = 1.31 
and a ’^corrected” version where the positive bulge produced by the 
local numerical disturbance (Figs. 6,13) has been removed. The dashed 
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line indicates . this ’region. As may be seen by comparison with Fig, 

21, the shape of this profile agrees qualitatively with experimental 
data. 

Figure 16 shows two skin friction profiles, one corresponding 
to a case approaching incipient separation, the other to that of a 
slightly separated region. It is observed that in the former case 
is much closer to x^^ than in the latter case, indicating a 
weaker interaction. The bulge that occured in the lower figure at 
the reattachment region (substituted by a dashed line) did not 
occur in the weaker case, indicating that the local numerical dis- 
turbance which apparently produces the bulges (Figs. 6,13,14), is 
associated only with those cases that involve a strong enough 
interaction to induce separation. 

Figure 17 shows the behavior of over the plate with con- 
stant suction velocity of approximately 5% of the normal velocity 
at the boundary layer edge. This suction is applied at the indicated 
region, starting at x^ and ending slightly downstream of the shock 
impingement point. lit is observed that all the numerical instabil- 
ities at the separated region have disappeared, resulting in a smooth 
profile. The separated region itself is smaller than for a similar 
case without suction (Fig. 14). The usual positive bulge at the 
reattachment region has also disappeared 

Figure 18 shows plots pf pressure and skin friction over 
the plate for the case = 1.31. It is observed that the region 
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of sharply decreasing coincides with the region of pressure 

rise downstream of The upper part of this rise together with 

the slight pressure plateau and the first part of the sharp rise 

in pressure occur over the separated region (C^ < 0) and further 

downstream both curves approach a constant downstream velue. 

Figure 19 is a plot of a typical streamwise distribution of 

boundary layer thickness over the plate, starting at x . There is a 

o 

sharp increase in thickness until a peak is reached at the shock 
impingement point, followed by a sharp decrease that forms the 
well known "neck", with a subsequent gradual increase, forming a flat- 
plate-like shape further downstream. (See also Fig, 12.) Note that 
at the shock impingement point the boundary layer is 50% thicker 
than at x^ . Note also that the curve is concave upward through 
X = 1.19 which explains why 6**' > 0 in the strong interaction 


region. 



VII. COMPARISON WITH EXPERIMENTS 

The experimental work of Hakkinen, Greber, Trilling and 
Abarbanel [4] has been used here for the purpose of comparison with 
the present numerical work. The experiments were conducted with an 
insulated wall, whereas the theory is for. a cooled wall. 

Figure 20 shows the result of a pressure profile comparison. 

It is observed that the pressure rise starting at is slightly 
sharper than the ope predicted by the present analysis. Further- 
more , the pressure plateau is better defined than the present work 
predicts. Also, the second rise in pressure downstream of the 
shock impingement point, while exhibiting the same behavior, is 
predicted somewhat upstream of the place indicated by experiment. 
Lastly, while the program predicts a rather abrupt rise to the 
downstream pressure, the experiment indicates a more gradual behavior. 
All in all, the comparison shows a reasonably good fit between 
theory and experiment despite the mismatch in wall temperature. 

Figure 21 compares skin friction profiles. It is seen that 
the decrease in towards separation is more gradual than predicted 
by the program, and, while the point of separation is almost the 
same, the present work predicts a smaller separated region than 
indicated by experiment. Both these predictions, i.e., higher 
with more rapid drop-off prior to separation and reduced separated 
region in comparison with experimental data, have been found also by 
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a 

■ 0 

Dwoyer [9], Dovmstream of reattachment the larg“e discrepancy 
between the present work and experiment is due to the large positive 
bulge that is caused by the local nuirferical disturbance at the first 

t 

reattached streamwise station*. In summary, the present work compares 
quite favorably with experiment through separation. 


9 


VIII. CONCLUSIONS 

The present work indicates that using the compressible boundary- 
|l layer equations coupled with an exterior Prandtl-Meyer flow does 

lead to an approximate solution to the problem of oblique shock 
waEve-laminar boundary layer interaction. 

The agreement with experiments, although not completely satis- 
factory, is nevertheless remarkable in light of the theoretical and 
numerical difficulties involved in the use of parabolic equations 
for the solution of an elliptic problem. 

It is noted that the application of a relatively small amount 
of suction in the shock impingement region significantly reduces the 
size of the backflow region. Introducing suction into the cal- 
culation also proved to be especially helpful in reducing the 
severity of the numerical instabilities. 

The use of a linear finite difference approach in the present 

work, while significantly simplifying the formulation of the equations 

and reducing computing time, might have nevertheless been part of 

the reason behind the local numerical instabilities that occured in 

critical locations in the strong interaction region. Employing a 

non-linear Crank-Nicholson method for the finite differences might 

have reduced or even eliminated these instabilities, and might have 

also permitted a reduction in the number of mesh points used in the 

strearawise direction. % 
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The numerical differentiation procedure for 6*", when performed 

at the first station downstream of shock impingement, is most probably 

responsible for the local numerical instabilities that occured at 

the region downstream of shock impingement. A practical way to 

eliminate this problemwould be to skip the iteration procedure on 

the (J + 1) station (in addition to J itself) and resume it at 
sh ^sh 

the (J + 2) station. 

X , » 

sh 

In summary, the present work does produce a gross picture for 
the behavior of the flow as a result of an oblique shock-laminar 
boundary layer interaction on a flat plate, especially for the 
pressure and skin friction distributions. However, the present work 
comes short of providing a reliable, detailed picture of the flow in 
the strong interaction region. 
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APPENDIX- A 

BASIC DEFINITIONS FOR QUANTITIES USED IN THE PROGRAM 
1. Displacement thickness: . r •: 



(Al) 


where is the x-velocity component obtained from an inviscid 


solution of the flow field and evaluated at the wall. Likewise for 


In the program a Simpson's Rule is used for the integration 
and the upper limit is determined by the simultaneous satisfaction 
of the, following two requirements: 


r 




■jUj 

•aij 

■j y 


< s 
<£ 


where e is a prescribed" tolerance. This determines the edge of the 
boundary layer for each x-wise station. 


2. Skin friction coefficient and the criterion for reversed flow: 



(A2) 
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Hence the condition 0, 

j 

follows : 




Y-o 


Ay 


indicating backflow, is determined as 


so that <0 


or 


is equivalent to 

<0 



0 



APPENDIX B 


NEWTON-RAPHSON PROCEDURE FOR THE PRANDTL-MEYER EQUATION 


The general Newton-Raphson equation is: 

V -y ffyj) 

^ur i 

In our case, 
and 


where 


d4> 

dx 



Now 


_ I r? _ r^fu vJl 
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(B3) 
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and hence 


d _M#o\afx“"od)r/C^““^ ) cC X 

dx*M^ Ml 

Order of magnitude analysis for the terms in Eqn. (B3) 

x=2.fi.c.,M,=:a): 
x/n^ - Q -T., 

X J(‘-/ ~ 

For a flat plate 

j££~ g.s~ 

t^tT*g 0.8g 

<i(-X yRe, 
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Also 
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and 
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Hence , 

Away from the immediate vicinity of x = 0 this term is negligible 
in comparison with the other two terms in eqn. (B3) . Hence 

dM 

00 

is negligible, noting that 


- 

Mo« 


J 1 

U’-S'* 

m 

dXcCM^ 


and the term in the brackets is typically of 0(1) for a flat plate 
So the Newton-Raphson procedure becomes 




APPENDIX C 


C.l Numerical integrations of 6*" 

Since the iterations on each j station begin only at a point 
which is one step downstream of the beginning of the strong inter- 
action, at that stage is usually larger than 4, typically 50 
(allowing the development of a weak interaction over 49 stations 
before triggering the strong interaction'' . Hence the integration of 
6*" employs an inversion of the f basic formula for the 

second derivative, using a back ^ scheme, namely 


d"-H 

i ~ 2. CM 

Based on this we have 



(Cl-1) 


These formulae, when applied for the first two times, use the in- 
compressible weak interaction values for 6* and 6*', namely 

^* = l.7eosJ^ t 5 ^*'= 


where 
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C.2 Second derivative of 6* 

c 

Special care must be exercised when numerically differentiating 
a variable, since differentiation tends to introduce instabilities. 
It was found that the best formula for 6*" is a simple backward 


scheme , namely 


i ” tA *r \*- 


(C2-1) 


C.3 Newton-Raphson Procedure for Improving (6*"). 

S I 

The basis for the convergence of a calculation at each x-wise 
station is the agreement between the guessed and calculated 6*" for 
that station to within a specified tolerance. Hence the Newton- 
Raphson procedure is 
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i-l 


(C3-1) 


where j stands for any x-wise station and c means a computed value. 
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C.4 Newton-Raphson procedure for improving (6*") 

S J . 


sh 


The basis for obtaining a successively better guess for 6*" 
is picking that branch which will result in the correct down- 
stream pressure with negligible pressure gradient. Hence the Newton- 
Raphson procedure is 



(s; 


1 

]t[ 


m 


T 


(C4-1) 

where is the inviscid downstream pressure (see Appendix D) and 

P[(5*" L] is the pressure at the downstream station, as calculated 
sh 

using the jth guess for 6*" . The starting values used in the pro- 

^sh 

gram are: 
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with = -0.1 


with K 2 = 1.5 



APPENDIX D 


FIRST GUESSES FOR 6*" and x 

J o 

^sh 

Since the downstream boundary condition demonstrates a nodal- 
type of behavior, extremely small changes in, (6*") 

sh 

produce large changes in the branching solutions. Hence 

the first guess for 6*" is important in order to start a branching 

. 

sh ■ ' , 

solution. A good initial guess enables the employment of small 
corrections to it as an effective control parameter. For incom- 
pressible flow over a flat- plate, 
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(Dl) 


from which 


e-’‘‘'=-aAsoel 




(D2) 


and 


In nondimens ional form this becomes 


(D3) 
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(D4) 
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APPENDIX E 


PROCEDURE FOR THE x ITERATION 

o 


The following method is employed for the successive improve- 


ment of X : 
o 


'Ck$^c 


(XoJl** 0(o)l * (Kz)c 


(El) 


(E2) 


where .< 1 to produce a smaller x^ and thereby a more compressive 
branch, and K 2 > 1 to produce a larger x^ and thus a less com- 
pressive branch. Once a correction has been done on x , the K's 

o 

are updated as follows: 




(K,)t^0.SS 


(E3) 




One of these new K^s is used in the next xj- iteration to produce a 

corrected either more expansive or more compressive branch, as 

necessary, until the right branch is found. This iteration is per- 

iteration, explained in 
"'x 

Appendix C4. 


formed in conjunction with the (6*")^ 
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Fig. 3 Effect of changing the location of the start of the strong 
interaction, x^, with fixed location of shock impingement. 
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flXIRL DISTANCE, 

Fig. 13 Typical static pressure profile. 
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flXlfll DI51flNC[. x/x^^ 


Fig. 16 Conparison between akin friction plots. Upper plot depicts 

case approaching Incipient separation. Lower plot depicts ^ 
case with snail separated region. <M •2.0, Re . • 3.89x10^, 
Pr • 0.70, • 1.0, dashed line not produced by* program.) 
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^^8* 17 Effect of suction on profile. 
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Fig. 20 Comparison between experiment and theory. Experimental 
points are from Hakkinen, Greber_, Trilling and Abarbanel 
(4J. (M^ » 2.0, Re^j^ - 4.02xlC)5, Pr - 0.72, p^p^ » 1.47, 

Sutherland Law, T * = 295®K, s* .» 110®K) . 

* . . SC ^ - 
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Fig. 21 Comparison between experiment and theory. Experimental 
points are from Hakkinen, Greber, Trilling and Abarbanel 
[4]. (M^ - 2.0, Re^^ - 4.02 x 105, Pr = 0.72, p^/p^ = 1.47, 

Sutherland Law, T * -= 295®K, = 110“K) . 
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